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Multiscale modeling methods for the analysis of metallic microstructures are discussed. 
Both molecular dynamics and the finite element method are used to analyze crack 
propagation and stress distribution in a nanoscale aluminum bicrystal model subjected to 
hydrostatic loading. Quantitative similarity is observed between the results from the two 
very different analysis methods. A bilinear traction-displacement relationship that may be 
embedded into cohesive zone finite elements is extracted from the nanoscale molecular 
dynamics results. 


I. Introduction 

C lassical fracture mechanics is based on the comparison of computed fracture parameters to their empirically 
determined critical values. Although this paradigm has been extremely successful for modeling crack growth at 
structural scales, it does not describe the fundamental processes that govern fracture. Ultimately, the understanding 
of events and processes that occur at length scales on the order of 10" 9 to 10" 3 mis needed to understand crack 
growth. 

Multiscale modeling provides an efficient means of interrogating deformation and fracture of metallic materials 
from the nicro- to the nano-scales. Many multiscale modeling strategies have been explored by mechanists in 
recent years; a few of them, coupling methods, quasicontinuum methods, equivalent continuum mechanics methods, 
and modeling with cohesive zone models, are discussed here. Coupling methods involve a priori dividing a 
problem into spatial regions based on the simulation technique to be used and tying those regions together with 
interface boundary conditions. For crack problems, a recently advocated approach combines ah initio quantum 
mechanics at the crack tip, molecular dynamics (MD) near the crack surface, and finite element (FE) continuum 
mechanics in the surrounding medium 1 " 3 

Quasicontinuum methods " employ a continuum framework with a finite element discretization wherein each 
element is formed by nodes that coincide with “representative” atoms . 5 In regions where the behavior of each atom 
is essential to understanding the problem, all atoms are considered as representative and coincide with nodes. In 
regions where a coarse mesh is appropriate, a finite element may encompass many atoms within its interior in 
addition to the atoms that coincide with the nodes. Also, it is assumed that the element constitutive relations are 
obtained directly from atomistic calculations. The Cauchy -Born hypothesis is used to relate homogeneous 
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deformations of the continuum to the superposed lattice to obtain the element strain energy by summing up the 
contributions of the interatomic potentials of the contained atoms . 4 ’ 5 

The equivalent continuum mechanics method 3,6 is similar to the quasicontinuum methods, but uses the meshless 
local Petrov-Galerkin (MLPG) method rather than finite elements. In general, meshless methods are developed to 
overcome some of the disadvantages of the finite element method, such as the need to interpolate discontinuous 
secondary variables across interelement boundaries and the need for remeshing in large deformation problems. 

In nultiscale modeling with cohesive zone models (CZM), 7 " 9 the CZM can be inserted within finite element 
models of microstructures, either within grains to model transgranular fracture or along grain boundaries to model 
intergranular fracture. Typically, the properties of the CZM are determined empirically or heuristically. However, 
an ideal CZM-based multiscale modeling strategy would use MD results as input to the CZM. Thus, the CZM may 
be used as a means to combine these two methods to bridge length and time scales and to provide the basis of a 
modeling approach that spans from nanometers to millimeters . 

In this paper, multiscale modeling is explored to study intergranular fracture in aluminum, with molecular 
dynamics and the finite element method (FEM) as the building blocks. First, MD and the FEM are discussed 
individually to highlight their potential in the proposed multiscale modeling strategy. Previous attempts to combine 
the two are also mentioned. Then, MD and FE simulations on a bicrystal model are performed to extract the 
constitutive behavior of the grain boundaries. Finally, a method of recasting this information in the form of 
cohesive zone models is discussed. 

II. Background on Elements of a Multiscale Modeling Strategy 

A. Molecular Dynamics 

To effectively model poly crystalline materials at the atomic scale, large numbers of atoms are needed. 
Approximations describing the potential energy of the interactions between neighboring atoms have been developed 
in the form of empirical and semi -empirical potentials. 10 These potentials are defined for atom-pair and many -body 
force interactions. An MD simulation is made tractable for large systems of atoms by treating each atom as a point 
mass, summing forces due to interactions with surrounding atoms, and using Newton’s Second Faw over a 
succession of time steps to obtain trajectories specifying atom positions and velocities. 

In the case of metals, metallic bonding is formed through delocalized electrons that are effectively shared 
between all atoms in the lattice. Thus, simple atom-pair potentials cannot be used to account for all of the mutual 
interactions between the electron clouds surrounding the atoms. Modifications, such as cluster potentials, pair 
functionals, and cluster functionals, have been made to pair potentials to develop general classes of interatomic 
potentials. The embedded-atom method (EAM), which uses a pair functional and treats each atom as being 
embedded in a field of electrons created by the surrounding atoms, has been successful in simulating the atomic 
bonding in metals. 11 

Considerable attention has been given to the atomic -level dynamic instabilities that occur during fracture. In 
homogeneous media, these instabilities may result in crack branching and roughness at the crack surface 
dramatically increasing the crack area and the energy dissipation at the crack tip. 12 In ductile poly crystals, such as 
face centered cubic (FCC) metals, the transgranular crack propagation taking place through the grains experiences a 
“dynamic brittle -to -ductile transition” as a result of the dynamic instability. The transition is associated with 
blunting of an initially brittle crack in an explosion of dislocations at the crack tip 12 together with crack branching 
along the available slip planes and surface roughening. 13 In intergranular fracture, in which cracks propagate along 
the grain boundary (GB) interface, the crack tip is channeled by the GB layer. This constrained condition may affect 
the crack dynamics and instability in the crack motion and yield significantly different behavior compared with 
transgranular fracture. 

B. Finite Element Modeling with Cohesive Zone Models 

Cohesive zone models approximate an interfacial traction-displacement relationship using bulk material 
properties and can model the appearance of fracture surfaces in a continuum. 14 Although the normal and shear 
components of the traction and displacement may be considered separately, as in reference % a more realistic 
approach may be to consider that the components do not act independently of each other. In the Tvergaard and 
Hutchinson 7 coupled cohesive zone model (CCZM), the normal and shear components of the traction and 
displacement are combined into single measures, x and X, respectively, so that the responses are coupled. 9 The 
CCZM given in reference 9 defines a traction potential, 15 
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where A. is a non-dimensional measure of the relative opening and sliding displacements (4 and §) as defined by 



8^ c and 8[ C are the critical values for the opening and sliding modes, respectively, and a fracture surface is assumed to 
have formed once X reaches a value of unity. The use of cohesive zone models in finite elements permits highly 
discontinuous displacement behavior to be exhibited in a single discretized model without remeshing while using an 
incremental piecewise linear solution algorithm. The displacement discontinuity may be defined in three separate 
regimes that depend on the relative CZM element nodal displacements. First, the CZM element functions as an 
elastic constraint to enforce node coincidence at the GB between connected continuum elements. Secondly, during 
strain softening of the opening of the CZM element, the effective stiffness is reduced as a function of the relative 
displacements of initially coincident nodes. Finally, after the boundary has opened and the effective stiffness 
becomes zero, further increments may cause local load redistribution resulting in separated surfaces coming back 
into contact. In this regime, the CZM element functions as a contact element to prevent interpenetration. 

As a consequence of the unique numerical features of the CZM element, a number of issues related to 
convergence exist that must be addressed prior to performing a successful analysis. One issue involves the 
magnitude of the initial “penalty stiffness” that is used to enforce the elastic constraint. This quantity is problem- 
specific and can lead to convergence difficulties if it causes the global stiffness matrix to become ill-conditioned. 16 
Also, if interpenetration is detected after a boundary has opened, the abrupt replacement of zero stiffness terms 
between surface elements with penalty stiffness values can cause divergence of the solution. Another issue pertains 
to the integration scheme used to compute the stiffness coefficients for the CZM element. In some instances, a 
Gaussian quadrature scheme can contribute to poor convergence compared to a Newton-Cotes scheme. 17 Solution 
convergence is also dependent on the size of the CZM elements compared to the extent of the process zone 
surrounding the fracture front. 18 An additional issue pertains to the possible extensive use of CZM elements within 
a model. If no a priori knowledge of crack paths is available, CZM elements may be placed between all continuum 
elements such that the crack path will be automatically determined by the predicted stress fields. However, if the 
initial or “penalty” stiffness is set to a small value, this can introduce an artificial flexibility into the model that can 
alter load transfer and displacements. 

C. Combining Molecular Dynamics and Finite Elements 

Attempts to extract relevant parameters for the decohesion law of a CZM element from atomistic (molecular- 
dynamics or molecular-static) simulations have been made in the last few years. 19 The approach in reference 19 is 
based on simulating the debonding of a flat interface with no initial flaw under a constant tensile strain rate 
perpendicular to the interface. In these models, the predicted mechanism for interface decohesion is not that of 
crack propagation, but rather of atom adhesion at the interface. Thus, the simulations reproduce the process of 
adhesion rather than that of mechanical fracture at the GB. Raynolds, et al. 20 use a similar setup for studying, by 
first principles calculations, the conditions and the mechanisms for initial void formation and crack initiation in 
adhesion curves of an NiAl-Cr interface. However, the resulting adhesion curves cannot be used directly to derive 
the constitutive decohesion laws for CZM elements because the boundary conditions at which a typical CZM 
element operates in a large scale FEM model are very different from the ones used in the MD simulations. 

III. Elements of a Multiscale Modeling Strategy 

The work presented in this paper is a potential step toward an integrated multiscale modeling strategy for 
understanding the connection between macroscopic fracture and the details of the underlying atomistic failure 
mechanisms. The range of relevant length scales span from the atomistic (10' 9 m) as shown in Fig. 1(a) to the 
macroscopic (>10~ 3 m) as shown in Fig. 1(f). Inherent to multiscale modeling across length scales is a trade-off 
between mechanistic detail and numerical efficiency. The mapping of the relevant physics between length scales 
necessarily involves a homogenization procedure to obtain an aggregate response over a representative volume of 
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Figure 1. Multiscale modeling with cohesive zone models. 

material. In general, this homogenization is performed at each increase in length scale and aims to preserve 
pertinent measures of fracture while reducing the number of independent variables. Within this modeling paradigm, 
various computational methodologies must be integrated to provide an effective homogenization. The particular 
homogenization used herein involves modeling the material behavior at the atomistic level using MD simulations 
and recasting the discrete interactions as local continuum field quantities through a cohesive zone modeling 
approach. 

The specific problem of intergranular fracture is considered wherein the MD simulations of a bicrystal (Fig. 1(a)) 
are used to obtain an atomistic response of crack propagation under applied loads. This response is then recast to 
obtain an averaged continuum traction-displacement relationship (Fig. 1(b)) to represent the cohesive interaction 
along a characteristic length of the GB (Fig. 1(c)). For the present stage of ongoing research, the resulting CZM is 
utilized as a bridge between atomistic and local continuum field modeling (Fig. l(d,e)). The determination of the 
dependence of the derived constitutive relation on length scale within a continuum representation is an ongoing 
research objective. 

A. Molecular Dynamics and Finite Element Analyses on Elastically Equivalent Systems 

A simple bicrystal model with a single grain boundary is used to study the behavior of a particular system as 
predicted by molecular dynamics and the finite element method. If these two fundamentally different analysis 
techniques can be shown to yield consistent results for this problem, then there exists a solid basis for the proposed 
multiscale modeling strategy. 

The MD simulation set-up used in the current analysis represents a thin bicrystal strip with periodic boundary 
conditions in all three directions. The atomic microstructure of the system is presented in Fig. 2. The normal to the 



4 

American Institute of Aeronautics and Astronautics 





L = 120 nm 

Figure 2. Initial atomistic snapshot of the MD system 
before crack propagation. 


strip (z-direction) for both crystals is along the 
[1 10] crystallographic axis. The two crystals, 

marked as “Crystal I” and “Crystal II”, are 
misoriented symmetrically at an angle of 89.42° 
relative to each other, forming a <1 1 0> £99 
symmetric tilt GB through which the crack 
propagates. The atomic structure of this high-angle 
GB in A1 is described in reference 21, and a complete 
discussion of crystallography is given in reference 22. 
The grain boundary energy, y G B, is estimated from the 
A1 potential to be 0.65+0.05 J/m 2 . 23 This large excess 
(i.e., above the perfect crystal) energy facilitates the 
GB decohesion. 24 The surface energy of the {5 5 7} 
GB plane, the plane along which the crack 
propagates, is estimated at y s = 0.945+0.005 J/rr? at a 
temperature of 0 K. The thickness (\|/) of the system 
(Fig. 2) in the z-direction is limited to 10{2 2 0} 
atomic layers (10/V2 lattice parameters), or 2.9 nm. 
This thickness is chosen because it is large enough 
(more than four times larger than the range of the interatomic potential 23 ) to prevent interference of the atoms with 
their periodic images along the z-direction, thus preserving the local three-dimensional physics in the system. The 
configuration allows the system size in the x- and y -directions to extend up to 120 and 97 nm, respectively, while 
limiting the number of simulated atoms to 1,994,000, thus allowing the simulation to be carried out on a modest 
Beowulf cluster. In addition, the <1 1 0> texture makes two of the four {111} slip planes in the FCC lattice 
available for dislocation glide and cross-slip events, as in a fully 3D space. 25 

The two narrower layers on both sides of the bicrystalline system (Absorbing Layer I and II in Fig. 2) serve as 
shock-wave absorbers 26 where a damping friction coefficient is applied to the atoms to absorb the phonon waves 26 
expected to be generated from the crack tips. The additional GBs formed between these layers and Crystals I and II 
act as absorbers for the dislocations that may be emitted from the crack tips during propagation. The inclusbn of 
these shock-wave and dislocation absorbers suppresses the creation of periodic images of all the disturbances 
emitted from the crack tips, thus avoiding the unrealistic influence these disturbances would otherwise have on the 
crack propagation. 

The simulated strip is prestrained by applying a constant external hydrostatic load (a** = a yv = c zz = a) at a 
constant temperature of 100 K. The temperature 100 K is chosen to suppress grain boundary and surface diffusion. 
After reaching mechanical equilibrium, the system size is fixed. The hydrostatic loading helps to eliminate 
plasticity effects not related to the crack, such as spontaneous dislocation nucleation from the GB, which could 
otherwise dominate the deformation in this (essentially) two-dimensional configuration. The crack is initiated by 
screening the interatomic potential between atoms on both sides of the Crystal I-Crystal II GB plane (Fig. 2) along a 
region of 5.7 nm along the middle of the GB. From this initial crack, the GB opens, and the crack starts to propagate 
in both directions along the GB interface. Figure 3 shows MD snapshots of cracks that propagate in a plate 
prestessed at initial hydrostatic loads of c = 3.5, 3.75, 4.0, and 4.25 GPa. Common neighbor analysis (CNA) 27 is 
used to identify atoms in different crystallographic states: FCC, HCP, and non-crystalline atoms (gray, red, and blue 
in Fig. 3, respectively). Atoms with more than 1/3 of their nearest neighbors missing are identified as surface atoms 
(green in Fig. 3). In each of the cases presented, the crack propagation is not symmetric in the -hr and -x directions 
along the GB (as defined in Fig. 2). The propagation in the -x direction produces two twin patterns (© in Fig. 3) 
that grow almost symmetrically in the two joined crystals and retard the -x direction propagation of the crack. Also 
observed in the snapshots are stacking faults (© in Fig. 3) and partial dislocations (CD in Fig. 3) that originate near 
the -x direction tip of the crack. By contrast, the propagation in the +x direction is more brittle, accompanied by 
only a few perfect dislocation emissions (© in Fig. 3). As a result, the crack propagates much further along the GB 
in the -hx direction. 

The corresponding FE configuration is presented in Fig. 4. The dimensions of the system are scaled to recover 
the proportions of the MD system. The model contains a built-in lenticular slit of varying relative length, 0.05 = l/W 
= 0.91, to study the evolution of the system at different stages of crack growth. Although the lenticular slit is not a 
crack in the true continuum mechanics sense, it is an analog to the MD configuration and hence is used to represent 
a crack-like discontinuity (hereafter referred to as a crack for simplicity). To avoid stress singularities in the FE 
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Figure 3. MD snapshots of cracks propagating in the MD 
system shown in Fig. 2 prestressed at initial hydrostatic loads 
of a = (a) 3.5, (b) 3.75, (c) 4.0, and (d) 4.25 GPa. 


system, the crack tips have a finite, but very 
small, initial radius that corresponds to an 
initial crack opening of approximately 1 nm, 
i.e., slightly larger than the range of the 
interatomic forces (0.67 nm). There is no 
elastic equivalent for the GBs in the FE 
system. It is assumed that the somewhat 
different elastic response of the GB layers 
does not alter significantly the crack behavior 
because of the relatively small volume ratio of 
the GB volume in the system. The matrix is 
assumed to be anisotropic linear elastic, with 
constants assigned according to the MD 
interatomic potential. Generalized plane strain 
is used, and the calculations are carried out 
under displacement control in the x- and y- 
directions. No dynamic or non-linear effects 
are incorporated, and the results are relevant 
for infinitely slow crack propagation under 
small strain. 

Plastic processes (twinning and dislocation 
emission 25 ) at the crack tips have a pronounced 
effect on the stress distribution around the 
growing crack. This effect is revealed by 
comparing the stress distributions obtained 
from the MD and the FE simulations of the 
two elastically equivalent models. The MD 
simulation necessarily incorporates all the plastic processes together with the elastic response of the bicrystal, while 
the FE model can be used to predict the elastic response separate from the plastic components of the solution. 

The MD simulation necessarily incorporates all the plastic processes together with the elastic response of the 
bicrystal, while the FE model can be used to predict the elastic response separate from the plastic components of the 
solution. The MD stress distribution is determined after the atomistic system reaches elasto-plastic equilibrium and 
the crack stops growing, while the FE stress distribution is determined for a static solution when the system is in 
elastic equilibrium. Two-dimensional (r-y) stress-maps for the MD results are created by averaging the local virial 
stress over a chosen volume of 6x6xlO/V2 lattice parameters in the x-, y-, and z-directions, respectively (2.43 x 2.43 

x 2.86 nm 3 volume including 864 atoms) and 
over 16 ps of simulation time. It should be 
noted that the \irial theorem, which is the 
basis for calculating the virial stress, provides 
the oldest and most frequently used expression 
for relating forces and motion within an 
atomic system to a continuum stress. In the 
limit of time and ensemble averages, the virial 
stress coincides with the Cauchy stress used in 
continuum mechanics , 28 

Figure 5 presents a comparison between 
the stress distributions obtained from the FE 
(Fig. 5(a-c)) and MD (Fig. 5(d-f)) simulations 
for the case of the small prestress of 3.5 GPa, 
which produces a 10 nm crack in the MD 
simulation at equilibrium (seen also in Fig. 
3(a)). The size of the crack is sufficiently 
small compared to the system size, and the 
elastic stress field calculated by the 

T-i* . t~i* i * . n j , ... . , corresponding FE simulation does not 

Figure 4. Finite element mesh ot model with built-in lenticular A , _ r t 

experience edge effects from the boundary 

crscK • 
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conditions and thus behaves as an infinite plate. The MD 
stress field captures the GB effect on the stress (© in Fig. 
5(d)). The MD results for stress distribution shown in Fig. 
5(d-f) for the +x crack tip are in very good quantitative 
agreement with the FE results shown in Fig. 5(a-c). This 
similarity implies that the crack propagation in the MD 
system in the +x direction is almost brittle with essentially 
no plasticity. For the stress distribution for the -x crack tip, 
there is a large difference between the MD results shown in 
Fig. 5(d-f) and the FE results shown in Fig. 5(a-c). This 
difference shows that the twinning observed at the -x crack 
tip in the MD system significantly alters the stress; the G xx 
and G xy stress components are largely relieved, while the 
o yy stress (® in Fig. 5(e)) is redistributed away from the 
crack tip along the twin boundaries (see Fig. 3(a)). This 
redistribution explains the very slow crack propagation in 
the -x direction of the MD model. As expected, the FE 
stress distribution is symmetric for the two crack tips 
because they are elastically equivalent in the FE model 
(Fig. 5(a-c)). 

For a crack tip propagating in an essentially elastic 
fashion (i.e., the +x crack tip in Fig. 5), consistent results 
are obtained by the MD and FE simulations, indicating that 
there is a solid basis for the proposed multiscale modeling 
strategy for brittle fracture modes. The extraction of a 
traction-displacement relationship for the +x crack tip from 
MD results to be used in larger-scale FE calculations is 
discussed next . 

B. Defining a Traction-Displacement Relationship from MD 

Grain-scale simulations that use CZM to study fracture (such as those presented in reference 9) typically use 
heuristically derived relationships and input values to define the CZM. In such an approach, these values tend to be 
only gross estimates of parameters such as GB yield stress; however, a cohesive formulation can be part of an 
effective physics -based approach if constitutive parameters are used that optimize the similitude between atomistic 
simulation and continuum finite element results. The model can be developed from results of the MD analyses 
allowing the physical insight of MD to be embedded in the more computationally efficient FE models. 

The traction-displacement function of a CZM describes how the traction x 
developed at the crack surface depends on the crack opening A, (see Fig. 6). In general, 
traction-displacement functions must be defined for both normal and tangential 
components of traction. Here, the special case of normal opening (mode I) under 
hydrostatic load is considered. The tensile component of the normal stress near the 
debonding GB interface from the MD simulation, (f yy , can be used as an equivalent of 
the normal traction x n in the CZM. To obtain the mode I traction-displacement 
relationship , the MD simulation box is divided into thin slices of length A, as shown in 
Fig. 7. A set of test volumes (the dotted areas in Fig. 7) of dimensions A x 2 A x \|/ 
centered at the debonding interface are defined. For each test volume along the GB 
interface, the normal opening of the crack 4 as a function of x is estimated, and the local normal stress G yy is 
averaged to obtain G s yy as a function of x. Thus, calculating <fyy(x) and 8^(x) for each test volume at position x along 
the debonding interface, one can develop the functional relationship & yy (8^), which represents the traction- 
displacement curve of a CZM volume element (CZVE) of size A x 2A x \\f. By definition, A must be small enough 
so that the stress variation over the CZVE is much smaller than the stress itself. In the present simulation, A= 1.9 
nm, which, with a thickness in the z-direction of 2.86 nm, gives a volume of 2 x 1.9 2 x 2.86 = 20.7 nm 3 and contains 
about 1250 atoms that contribute to the stress-displacement response of one CZVE. There are 63 CZVE along the 
120 nm long GB interface in the system (Fig. 2). During the simulation, the c fyy(S^) state is scanned every 4 ps for a 
period of 200 ps, resulting in 3000 values from which to determine the (fyy(S^) dependence for each simulation run. 



Figure 6. A general 

traction-displacement 

relationship 
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Figure 5. Stress contours from the FEM model and 
the corresponding stress maps from the MD model. 
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Figure 7. Slicing of the system volume in the MD 
simulation for extracting the parameters for the 
cohesive zone volume elements (not to scale). 



& n , imi 

Figure 8. Surface stress vs. crack opening curves, 
<3 S jj(5„), characterizing the propagation of the cleavage 
tip in the +x direction for three prestrain loads. 


Fifteen hundred 0^(2^) values are taken from the CZVE on the right half of the model for the case of 4.25 GPa 
prestress and are presented as points in Fig. 8. The spread of the data points shows that at this small length scale of 
2 nm, there is a large dispersion around the mean traction-displacement curve. These points represent the behavior 
of the crack as it propagates in the +x direction (see Fig. 5). The +x direction tip propagates more than 50 nm (Fig. 
3), causing 25 CZVE to experience a complete transition from a fully closed to a fully opened state. The c ^(2^) 
curve for the 4.25 GPa prestress is extracted by taking a moving average of the <3^(2^) values over 0.0 < <4.0 

nm. This process is repeated for the cases of 3.75 GPa and 4.0 GPa prestress, and the resulting <3^(2^) curves are 
also shown in Fig. 8 (note that the corresponding <3^(2^) values are omitted from the figure). The curves in Fig. 8 
show a bilinear type of constitutive relation for the CZVE similar to the CZM suggested by Camacho and Oritz 8 and 
show that for the case of crack propagation in the +x direction for a Z99 STGB in FCC metals, this bilinear model is 
a good assumption. In addition, Fig. 8 shows that <3^(2^) is a weak function of prestress of the system over the 
range 3.75 GPa <o< 4.25 GPa and is approximately c p = 5.0 GPa at 8^ =0.4 nm. The full opening 

(debonding) of the interface happens at d c n =2.5 nm, when & yy becomes zero. 


IV. Concluding Remarks 

The present multiscale modeling strategy for the analysis of aluminum micro structures is based on embedding 
atomic -scale information into finite element simulations by means of cohesive zone models. Analysis of crack 
propagation in a simple bicrystal model by means of both molecular dynamics and elastically equivalent finite 
element analysis reveals a quantitative similitude in stress distribution for cracks propagating along a grain 
boundary. By examining the tensile component of the stress and the crack opening from the atomistic results, a 
bilinear traction-displacement relationship is developed that incorporates the behavior of approximately 1250 atoms 
into one equivalent cohesive zone volume element. 
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